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SYNOPSIS 

DETEMBTATION OF TEMPERATURE IROFILES IN QUENCHING- BY NUMERICAL 
METHODS AND VERIFICATION BY EXPERIMENTS 

H.K. Kanodla 

Master of Technology 
Department of Mechanical Engineering 
Indian Institute of Technology , Kanpur 

The present work is concerned with the prediction of temperature 
profile during the quenching of a material- Hie surface heat transfer^ 
coefficient during quenching vary s ignif ic antly in the film boiling and 
nucleate boiling regim* and also depends strongly on the temperature 
difference. Hence the problem for one dimensional heat conduction 
reduces to the solution of a partial differential equation in two 
variables with a nonlinear boundary condition. This is converted to 
nonlinear singular Volterra integral equation. Hie integral equation 
is solved numerically by two methods - the successive iteration and 
separable kernel method. The numerical solutions are compared with 
experimental results and results available in the literature.' It is 
shown that the separable kernel method is more accurate if sufficiently 
large number of terms are considered in the expansion of the singular 


kernel 



CEAFTER 1 


DWRODUCTiai 

Heat treatment is an important method used in industries to 
increase the quality, reliability and durability of engineering 
components by imparting the desired mechanical properties. The desired 
mechanical properties can be increased hardness, wear resistance, 
ductility, toughness etc. Heat treatment processes can be broadly 
classified into the following categories: 

1. Thermal Treatments: Where there is a change in type and 
proportions of phases present in the solid or concentration and 
distribution of crystal defects. 

2. Therm ochemical Treatments : 'Where there is a change in chemical 
composition of surface layers (e.g. carburising and nitriding). 

The thermal treatments can be subdivided into four sections 
depending upon the structure resulting from the treatments. They are: 

a) Annealing: Heating to and holding at a suitable temperature 
and then cooling at a suitable rate, for such purposes as reducing hard 
ness, improving machinability , facilitating cold working etc. 

b) Quenching: Hardening a ferrous alley by heating it above 
the phase transformation temperature and then cooling rapidly so that 
some or all of the austenite transforms to martensite. 

c) Tempering: Reheating a quenched hardened alley to a 
temperature below the transf ormation range and then cooling at any 
desired rate. 
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d) Aging ? A change in. the properties of certain metals that 
occurs at ambient or moderately elevated temperatures. 

The main aim of quenching and tempering is to increase the 
hardness and wear resistance, retaining a sufficient toughness at the 
same time. Quench hardening is discussed in detail hereafter. 

As mentioned above, quench hardening is one of the thermal 
heat treatment processes, where concentration and distribution of 
crystal defects play an important role. The mechanical properties 
of a crystal are determined by the number and mobility of disloca- 
tions contained in them. Point defects exist in crystals in thermal 
equilibrium and may contribute to mechanical properties through 
diffusion, however, remarkable effects are caused by nonequilibrium 
concentration of point defects. If this nonequilibrium concentration 
of point defects is achieved by rapid cooling from high temperatures, 
the resulting hardening is called "quench hardening". 

Hardening in the ferrous alleys is a result of supercooled 
transformation of phase. This can he explained by cooling curves on 
TTT (tine, temperature and transformation) diagram, Fig.1 . The TTT 
diagram is a curve, which gives the percentage and structure of decom- 
posed austenite corresponding to the temperature and time. If 
canponents are cooled rapidly, austenite does not have enough time 
to decompose at a high temperature with the formation of ferrite- 
cementite mixture. It is supercooled and is transformed into martensite. 
Martensite is a hardened steel structure. 
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She most extensively used method for hardening is quenching 
in a single medium. The main disadvantage of this method is crack 
formation due to simultaneous influence of transformation stresses 
and thermal stresses (caused by high rates of cooling in the temper- 
ature range of martensite formation). Hence, sometimes two quenching 
media are used. The second quenching media is less intensive. When the 
component is transferred from one quenching media to another quenching 
media, it may or may not he held at a constant temperature for s cme- 
time (Pig.2). 

For a given design of a component, factors life nature and number 
of quenching media, initial temperature of the component, temperature 
of quenching media and holding time in quenching media, depend upon the 
material chosen and desired mechanical properties. All this requires 
an extensive knowledge of the rate of cooling of both the surface 
and interior of the component during quenching. 

In quenching , heat flows from component to the quenching media 
at the interface, known as surface. Heat transfer at the surface 
mainly occurs ty convection and/or by radiation. For unsteady convective 
heat transfer, boundary conditions of the fourth kind is used i.e., the 
interface temperature of both the media and component are same and 
boundary layer theory is applied. If the quenching media :.a large 
its temperature can be considered uniform and heat flux is proportional 
to the excess -temperature C— wer fluid temper' ture ). Hie propor- 

tionality constant is known as surface heat transfer coefficient which 
may be a function of excess temperature. 




(a) Quenchrng m single (b) Stepped quenching 
medium ■ 



medium without 
time lag 


A c ^ Recry stall sat ion temperature ■ 
M - Martensite formation -STARTS 
Mj - Martensite formation - STOPS 

Methods 'of quenching 



Convection I Um tee to . -bos ling \lmnmsm FEm boiling 



Typical cooling curve for -a pool 
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As shown in the plot (Fig. 3), heat transfer coefficient is a 
strong function of excess temperature of surface. The entire region 
may he divided into four distinct stages of different intensity: 

a) At the first stage (legion A) a thin vapour film surrounds 
the hot metal. Cooling proceeds by film boiling. Cooling rate is 
slow and is determined mainly by radiation. 

b) In this region (Region B) cooling proceeds by transition 
boiling, a stage in between the film boiling and nucleate boiling. 

c) Hie vapour film breaks up and liquid boils with bubble 
formation on the surface of the metal being cooled. Since all cooling 
is accomplished by vapour generation, this is the fastest stage of 
cooling, known as nucleate boiling (Region c). 

d) Below saturation temperature, cooling is much slower as 
heat is extracted mainly by convection (Region D), 

The earliest analytical work for transient cooling was done 
by Fourier (1878) with simplest boundary conditions, namely of constant 
surface heat transfer coefficient. Results of this work have been 
tabulated and plotted for practical use by Curney-Iurie (1923) and 
Heisler (1947)* Though these results are presented for one dimensional 
heat flow only, these can be extended to three dimensional heat flow 
also. 

In hardening, cooling proceeds through all stages of heat 
removal. Hence surface heat transfer coefficient charges rapidly as 
surface temperature goes down. Hence, analysis given by Fourier ( 1878 ) 
cannot be applied directly. In order to use Heisler charts we need to 
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assume seme average constant heat transfer coefficient. However, 
there is no guide-line for the choice of this average heat transfer 
coefficient. 

For obtaining detailed history of temperature at ary time and 
coordinate during quenching, one has to go for mathematical analysis. 

The problem of transient heat conduction in a solid with surface heat 
flux as a function of temperature is a nonlinear problem. The nonlinear 
transient heat conduction problem can be reduced to the solution of a 
nonlinear singular Vol terra integral equation. Numerical methods are 
used to overcome the nonlinearity and singularity. 

The present study is concerned with a general nonlinear 
transient heat conduction problem resulting from a nonlinear surface 
flux. In general, solid is subjected to surface heat flux varying 
with various stages of boiling and natural convection. 

Ihile the literature dealing with the particular cases of 
convection, radiation is extensive, the unsteady problem of heat 
conduction with varying heat flux at surface due to various boiling 
stages has received little attention. Transient heat conduction in a 
plate subjected to pure thermal radiation on one face and combined 
forced convection and radiation on the other has been studied by 
Birka (19 66). The problem was reduced to the solution of two simul- 
taneous singular nonlinear Volterra integral equations. The two 
equations were replaced by a system of nonlinear algebraic equations and 
were solved for one set of parameters over a limited time range. 
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Crosbie and Viskanta (1966) presented an iterative numerical 
technique for solving nonlinear Volterra integral equations. The same 
method was applied to solve the problem of transient heating or cooling 
of a plate by combined convection and radiation ( 1968a) . They compared 
their results for one set of parameters with the results of Vidin and 
Ivanov (1965) and results obtained by finite difference method. 

Vidin and Ivanov (1965) transformed the linear energy equation 
and nonlinear boundary conditions into a nonlinear partial differential 
equation with a linear boundary condition. The resulting nonlinear 
partial differential equation was solved approximately by iteration. 

A simplified method which approximates the kernel of integral 
equation by a separable kernel has also been presented by Crosbie and 
Viskanta (1968b). This method was improved and applied to transient 
cooling of sphere in oil by Crosbie and Banerjee ( 1 97 3 > 1 974 ) • 

The present work is concerned with the development of a 
computer programme for obtaining the history of temperature vs. time 
plot for a quenched component. The computer may be regarded as a 
simulator and is used, in the same manner as an experimental rig of 
great flexibility in which the effect of changes in different parameters 
may be investigated. A quenching process can be finalised with some 
certainty before we go for actual testing, with their problems of 
high cost and time delay. 

Some experiments have also been performed and results compared 
with numerical solutions. Although, because of laboratory limitations, 
experiments for very high initial body temperature could not be performed. 
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The present work can be extended for different positions of 
components. Also, one may consider the variation of heat transfer 
coefficient with the coordinates. Though, this will be highly nonlinear 
problem. Complicated shapes, other than infinite plate, infinite 
cylinder and sphere, have to be solved by finite element method or 
finite difference method. 



CHAPTER 2 


ANALYTICAL FORMULATION 


in the present work only one dimensional geometries are consi- 
dered. Physical and mathematical models for infinite plate and infinite 
cylinder have been considered and their solutions discussed here. 


2.1 INFINITE HATE 


2.1.1 PHYSICAL MODEL s 


Consider a slab of thickness 2L with uniform initial temperature » 
Tf* Ik® slab is suddenly quenched in a cooling medium (also known as 
quenchent) at temperature T f , such that T f < db* Hie following 
assumptions are made in the analysis? 

1. The conduction heat transfer is one dimensional. 

2. The plate is isotropic> homogeneous and opaque to thermal 
radiation. 

5. The physical properties are independent of temperature. 

4 . Quenching media temperature is not a function of time . 
and 5. Hie plate does not contain any heat sources or sinks. 


2.1 .2 MATHEMATICAL MODEL? 


The governing 

3 T(x»t) 
3 t 


differential equation 
= a V 2 T(x,t) 


reduces to 


* 


9 T(x,t) = 3 2 T(xtt) 

3 t 3 x 2 


( 2 . 1 ) 


( 2 . 2 ) 
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for one dimensional heat flow. This partial differential equation has 
two boundary conditions and an initial condition. As problem is symmetric 
about the center plane, which is the origin of the coordinate system, 
the boundary conditions at x = 0 and x = L are 

2-h°Ul . 0 (2. 3 a) 

3 x 

a T(L,t) = _ i > 

3 x K 

and initial condition is 
T(x,0) = I £ 


(2.3b) 


(2.4) 


First boundary condition (2.3a), specifies that the temperature 
gradient at the center is zero, while second boundary condition 
(2.3b), is a statement of energy balance at the surface. To simplify 
the analysis, eqn.(2.2) is non-dimensionalised in terms of dimensionless 
variables and then the eqn.(2.2), reduces to 

38 (x, Fo) = 3 6(x, Fo) \ 

3 Fo 1 ? (2.5 ) 

3 x 


while the boundary conditions and initial condition reduce to 

= 0 (2.6a) 

q { 6 (1, Fo)} L 


38 (0, Fo) 
3 x 

ae(i, fq) 

3 x 


K 


(2.6b) 


and 9 (x, 0) = 1 


(2.7) 


I 
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•where 


T - T 


f 

0 - V % 

(2.8a) 

- _ X 

X " L 

(2.8b) 

Fo - “1 
x2 

(2.8c) 


2.1.3 SOLUTION TO THE PART IAL DIFFERENT IAL EQUATION 

Since the second boundary condition (2.6h) is nonlinear, none of 
the conventional methods are useful for solving the above partial 
differential equation (2.5). This equation can be reduced to a non- 
linear, singular Yolterra integral equation by different methods, as 
shown by Luikov (1968). In the present work this equation is solved 
by Laplace Transform method. 


Taking the Laplace Transform with respect to the variable Fo, 
partial differential equation (2.5), reduces to a total differential 
equation 


_ <i^F ( Xj s ) 

s 0(x, s) - 0 (x, 0) = 

d x 


where 

00 

F(x, s) = J exp(-s Fo) 6 (x, Fo) d Fo 


Boundary conditions (2.6), can be written as follows; 

d 9~ (0, s) = Q 
dx 

, dF (1, s ) F^(i, s)>'L 

and — L 

dx 


(2.9) 


( 2 . 10 ) 


(2.11a) 


(2.11b) 


K 



n 


Substituting the value of q (x, 0) in total differential equation. 
(2.9 )> from initial condition (2.7) and rearranging the equation, we 
get 

a 2 ? (*, al . „ 9 ( ;, s) ( 2.i2) 

dx 

General solution of this equation is 


9 (x, s) = A exp(/ s x) + B exp(-/ s x) + 


(2.13) 


where A and B are arbitrary integration constants. Using boundary 
conditions (2.11 ), for finding the values of A and B, general solution 
can be written as 


0(x, s) = - 


q{0 ( 1 » s ) } L exp(/s x)+exp(- /s x) ^ 


K 


yfe exp( /s) - exp(- /s) 


+ — 
s 


or 


_ q{^" (l , s ) } L cosh x / s 

e(x, s ) n 


K 


/ s sinh/ s 


(2.14) 


(2.15) 


Using the relation from leplace inversion tables given by Doetsch ( 1961 ), 


-1 r cosh x/s "*] 


/s sinh / s 


2 2 

1+2 l exp(-n Fo) cos { 2n tt( 
n=1 


x + 1 


)> 


(2.16) 

and using convolution theorem, general solution of differential equation 
( 2 . 15 ) can be written as 

x Fo 00 ? 2 

6 (x, Fo) = 1 - ?? J [1+2 £ exp{~n it (Fo -t)}cos{ mr(x+l)}] 


n-1 


(2.17) 


q {0(1 » x )> & t 
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or 

- 1 Fo 00 n _ 0 9 

0 (x, Fo) = 1 - J [1+2 J (-1; cos(nit x)exp( -n tt (Fo-t )} ] 

h 0 (1 ,T ))© (1,T ) d T (2.18) 

where h is heat transfer coefficient, it may he -written as 

<1(0(1, T )> - h{0 (1,T )>0(1,T ) = h ref 50 (1,T ) }0 (1, x) 

(2.19) 

where h re £ is a reference heat transfer coefficient. It may he an 
initial heat transfer coefficient or a maximum heat transfer coeffi- 
cient, whichever is convenient, h is a nondimens ional parameter. 
Substituting h hy h h in eqn.(2.18), we get 

Fo “ n - 2 

0 (x, Fo) = 1 - Bi / [ 1+2 l (-1) cos(nux) exp { -n it 2 (Fo-t )} ] 

r 1 o n=1 

h (0(1, r)}0(l,e )d T (2.20) 

•where is Biot no. corresponding to reference heat transfer 

coefficient. This integral equation is a nonlinear, singular Volterra 
integral equation and for surface this equation reduces to 

0(1, Fo) = 1 - Bi ref / [1+2 £ exp (-n 2 tt 2 (Fo^ t)} ] h{e(l, t) } 

o n=1 

0 (i , x ) dt (2.21 ) 

The kernel within the integral sign of eqn.(2.20) is known as Hheta 

function 0| (Appendix i). 

00 

0 . ( -H3 1+2 I (-l) n cos(nTTx) exp ( (n 2 'tt 2 (Fo-t )} (2.22) 

42 n=1 
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2.2 morns cylikier 

2 .2 . 1 PHYS IC A! MODEL s 

Consider a cylinder of diameter 2E and uniform initial temper- 
ature The cylinder is suddenly dipped in a quenchent at temper- 
ature such that < T^» Same assumptions are made as in the case 

of slab (article 2*1.1 ) and similar analysis is done for cylinder. 

2.2.2 MATHEMATICAL MODEL* 

The governing differential heat flour equation 

$) . -a V 2 T ( r ,t) (2.23) 

3t 

reduces to a nondimens ional form 

98 ( r, Fo) = 3 2 9 (x , Fo) [ 1 30 (t , Fo) ( 2 . 24 ) 

3 Fo 3 r 2 r 3 r 

for one dimensional heat conduction. Boundary conditions are 

39 (0, Fo) _ o 
3 r 

59 (1, Fo) = q{ 9 (1, Fo) 1 H 

a; “ E 

and initial condition is 
9 (r , 0) = 1 

•sdiere r = r/R 

This eqn. (2.24) with nonlinear boundary condition (2.25b) is also 
solved by Iuikov (1968), but here it is solved by using laplace 
Transform method. 


(2.25a) 

(2.25b) 

( 2 . 26 ) 

(2.27) 
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Taking Laplace transform of eqn. (2.24) ?ath respect to the 
variable Fo and using initial condition (2.26), partial differential 
equation reduces to total differential equation 

si + i s Q( r ? s) - -1 (2.28) 

d r2 r dr 

The general solution to the eqn. (2.28) is vrritten as 

9 ( r, s) = A I Q (r /s) + B K Q (r/ s) + (2.29) 

where A and B are arbitrary constants of integration and I and are 
modified Bessels functions of 1st and 2nd kind respectively and of 
zeroth order. As 


r 0, K Q (/s r) *+ co (2.3O) 

"While physically, temperature of the component cannot be infinite, 
therefore, constant B should be zero and solution reduces to 


0(r ,s) = — + A I Q (r /s) 


(2.31) 


Value of A can be calculated frcm boundary condition (2.25b) 
Theref ore , 

K / s Ij ( /s) 


9(r, »)•'•*■• 


(2.32) 


Using Appendix II and convolution theorem 

Fo 


" J o( p n 


e(r,Fo) » 1 - | / [ 2+2 I — exp { - p£(F 0 ~T ) >]q{6(l, t) }dx 

o n=1 o^ 

( 2 - 33 ) 


or 


Fo « J 0 (p n r) 


e(r,Fo) o 1- Bi^- / [2+2 £ rw — t— exjf -P 2 n (Fo-T )}]. h{8(l , T )} 

ei o n=l o v n' 

8 (1 , T )d T ( 2 . 34 ) 
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where y n are roots of characteristic equation 

J« 5 C y ) = Ji( u) = o (2.35) 

Solution eqn.(2o4) also is a nonlinear, singular Yolterra integral 
equation and for surface reduces to 

Fo 00 

0 (l,po) = 1- Bi f [2+2 l exp{-y 2 (Fo- t ) }] 

o n=l 

MO (1,t ) }0 (1,T )d t ( 2 . 36 ) 

Sqns. (2.2l) and ( 2 . 36 ) can be written in a ccmbined form 

Fo co 

8 (1, Fo) = 1-Si / [ $ + I B exp (-y 2 (Fo-t ) >] 

ei 0 n=1 

h{e(l,T )} 0 ( 1 , x)dT (2.37) 

where and 8 n are constants and y n are eigenvalues. Values for 
f^, 8 and y are given in Table 1 for three simple ge one tries. 

Nonlinearity of eqn.(2.57) does not allow any closed form 
solution. Hence, we resort to numerical methods to solve this equation. 
Numerical methods have been discussed in the next chapter. 
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Tab le 1 


S.ITo. 

Geometry 



6 k 

X k 

M o 

*1 

i. 

Infinite 

Plate 

1 

2 

k 

1/3 

1/45 

2. 

Infinite 

Cylinder 

2 

2 

J-j ( ^ jj) = ® 

1/4 

l/96 

5- 

Sphere 


5 

2 

tan ^ k 

1/5 

1/175 





CHAPTER 3 


NUMERICAL TECHNIQUES 

It was observed in the last chapter that the Volterra integral 
eqns. (2.20) and (2.34) were singular and nonlinear. Two different numeri- 
cal methods for solving equations of this type have been suggested by 
Crosbie and Viskanta (1968a) and Crosbie and fenerjee (1973)* These 
methods are discussed here. 

3.1 SUCCESSIVE APPROX 3MATICN METHOD 

This method was suggested by Crosbie and Viskanta (1968a). The 
nonlinearity in eqn. (2.20) arises due to nonlinearity of the boundary 
condition (2.6b) at the surface. For the surface eqn.(2.2C>), reduces 
to 

F° co ? 

0 (l,Po) = 1-B1 {[1+2 T exp{_ n it (Fo-x )}1 h{ e(l> T ) >6(1 »T )dx 

rei 0 n=1 

(3*1) 

Once eqn. (3*1 ) is solved for a dimensionless tine T, such that 
0 < x <Fo, then eqn.(2.20) can be solved for any coordinate x and time 
T , because 9 (l,T ) and hence each term on right hand side of eqn.(2.20) 
will be known. The above equation can be solved by successive approximation 
method. An initial approximation is made for 8 (l, t) and then it is 
improved by successive iteration, (k+1 )th approximation can be given 
as 

Fo 00 2 2 

8( 1 ,F°)k +1 = 1 “ ;Bi ref / [1+2 ^ exp{ -n w (Fo-t )> ] 

* Sted.T^e (i,t ) k d t 


(3*2) 
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For numerical solution, tine period 0 <x < Fo is divided into 
m small intervals, such that 


0 <*-! < t 2 < t 5 < ••*< V! < F ° 


(3-3) 


Then for the first interval 0< t < t^ , eqn.(3*2) can be re-written as 

1 °° 2 2 

0(1, t 1 ) k+1 = 1- Bi ref J [ 1+2 exp{-n it (t.,~x)}] 


h{6(l, T^} 0(1, T ^ d T 


(34) 


This equation can be solved directly by successive approximation 
method. For second interval, 0< T< ± t integral of eqn.(3*2), can 


be broken into two parts 


0 


(l,t 2 k +1 - 1-Bi ref I [1+2 I exp { -n 2 7r 2 (t 2 -T )}] 

0 n=1 

fi" {6(1, x )> 0 (1,T )d T 
^2 00 

- Bi__„ I [1+2 I exp {-n 2 F 2 (t 9 - t)>] 


ref 


n=1 


S {0(1, T ) k } 0 (1,T di 


(3-5) 


For first integral of equation (3*5 )> function 0 (l,T ) is known over 
a period of 0 < T< t^ from previous calculations. Hence eqn.(3*5), 
can be solved by assuming an initial approximation for 0 (l, t) in the 
time interval t] < x < t 2 and improving it by iteration. Similarly 
(k+1 )th approximation for mth interval is given by 

*m-1 00 2 2 

0(1 ,F°). ^ = 1-Bi . f [1+2 7 exp { -n it (Fo-x )>] h{0 (1, x )>0 (l,x)d t 

* +1 ref 0 n=1 


Fo 


-Bi r f / [ 1+2 I exp { -n tr (Fo- t )}]h {0 (l , x ^ 6 (l ,x ^ dx 
Vl ^ (3.6) 
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She choice for initial approximation is of great importance. 
Convergence -will he faster if initial approximation is closer to exact 
solution. This will save computation time. Crosbie and Viskanta (1968a) 


used the solution of the following equation as an initial approximation. 


(l,Fo) = 1-Bi j [ I42' l exp { -n ir 2 (po- t ))] 
rel 0 n=1 


h {6(1, T ))e (l,r )d T 


- Bi h{6(i_p 0 )} 0(1, Fo) 


/ [ 1+2 l exp {-n 2 Tr 2 (Fo- T )}] d T 
. 1 *-1 

m (3.7) 


The abo-ve equation can be solved for9 (l,Fo) as terms contained 
in the both integrals are only known terms and hence integrals can be 
evaluated. It can be seen that eqn.(3»7)> represents a first order 
approximation to eqn.(3.6), by neglecting the variation of h{0(l,T )}0(l,T ) 
in the time interval t m _i to Fo. On the other hand, when time interval 
tuj.i to Fo is small, value of 9(l,t^_i) can be used as an initial 
approximation for 9 (l»Fo) * This is being used in the present work; 

The kernels of integrals in eqns. (3»4)> (3*5) and (3*6) are singular 
because infinite series does not converge at t= Fo and summation of 
infinite series tends to infinity. This singularity of kernel can be 
removed by Poission summation formula (see Appendix III) 


I4-2 £ exp{ -n tt (Fo-t ^ = [ 1+2 J exp { -n/( Fo-x )}]{7r (Fo-t )} 

n=1 n=1 

(3.8) 

and after substituting it in eqn.(3*6), integral is evaluated by "the USe 
of modified Gaussian integration formula, Bevis and Polonsky ( 1965), which 


states 
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J f(x )(Fo- t) s dr = 2(Fo-t 2-1 ) z l f( Ti) (3.9) 

m-1 i=1 

where 

T • = t i + (Fo- t 1 )(l-X. 2 ) 

1 m-i v m-1 v x ' 

and = Gaussian weightage of order (2N). 

Substituting eqn.(3*9), in eqn.(3«6), we get 

^m-1 os 

9 (l,Fo) = 1-Bi ~ / f 1+2 l exp{ -niT (Fo-t )>] 

rei o n=1 

E{0(1,t)}0 (1,t )d t 

“ 2Bi ( A t/ir F l ,^h 1+2 l exp{ -n 2 / (At x? ) > ] 
rei i=1 n=1 1 

h{e (i,x i )>e (1, x ± ) (3.10) 

where A t = Fo-t .. . 

m— 1 

3.2 SEPARABLE KERNEL METHOD 

A method which is simpler than the successive approximation 
method was also suggested by Crosbie and Viskanta (1968b). This method 
is known as the 'Separable Kernel Method' . The kernel in the integral 
eqn.(2.37) is represented by an infinite series. 3h the separable 
kernel method, the infinite series is truncated at a suitable point, i.e. 
symbolically, the kernel of eqn. (2.37) 

00 2 

f = 3 0 + 3 n exp{ - (Fo-t ) } (3-11) 

is approximated by a separable kernel 
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fa=B + f B n exp{ -p*(Po-T ) } ^* 12) 

13=1 

Later this kernel was modified by Crosbie and Banerjee (1973), for 
taking into account the other terms of infinite series. Modified 
kernel is 


M-1 


fa = 


B 0 + I B n exp ( - w n (Fo- t)}+ e 8 exp { -b^Fo-T ) } 
n=l 

(3,13) 


Then the nonlinear singular Volterra integral eqn.(2.j57)> can be 
written as 

Fo 

0 (1 ,Fo) = 1 - Ei ref f 0 ^ h {6(i, T)>0 (i, T )d T 

M-1 Po 2 

“ I Bi ref B n q / exp( -P n (?o-T ) > 

h{ 6 (1 ,t )> 0(1, T )dx 
Po 

- Bi ref e M 6 M / exp { - b . (Po-x )}h{0 (l,x )} 6(1? )dx 

(3.14) 

and solution to the above equation can be split into (M+1 ) parts 

0 (l,Po) = 0 O (1, Po) + 0 .,(1, Po) + ... +0 M (l,Po) (3-15) 

where 

Po 

6 0 ( 1 » F °) ■ 1_ ^ref e 0 f E{0(l,x)}0 (l,x )dx (3.16a) 

0 

Po 2 

9.,(l,Po) = -Bi f 6 1 J expf-V-j (Fo- T ) > 

rex P 


h {0(1, T)>0 (1,T )d X 


(3.16b) 
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Fo 

® M-1 ^ 1 ’ F °^ = " Bl ref B M-1 ^ exp { - u M f 1 (Fo- t ) } 

Mb 0 » x)} 6 (1, x)d T 

and (M+1 )th term can be written as 

Fo 

9 » Fo ) = " Bi ref £ M B M I expf ' b M (p °- T > } 

0 

S{0(1,T )} 6(1, T )d T 


(3 . 16c ) 


(3.16a) 


Differentiating above equations with respect to Fo once* we will 
get (M+1 ) simultaneous differential equations connected by eqn.(5»15)* 
The set of differential equation is given as 


d 6 (l »Fo) 

= 1- 8 0 Bi ref h{ 0(l,Fo)>B(l,Fo) (5.17a) 

&Fo 


d6 1 (l,Fo) 
dFo 


- - 6 1 Bi rgf h {6(l,Fo)}0 (l» %)- V ^ 0^(1 ,Fo) 


(5.17b) 


d0 M-1 (1 ’ F °) - 2 

55 ' ' 6 M-1 ^ref h{6 Cl > F0)} 8(1 ’ P0) ’"h Vl 0 ’ 1 ^ 

(5.17c) 

and 

dB M (l, Fo) 


dFo “ht 6 * Bi ref h{0(l,Po)}e (l,Fo) - (1,Po) 


(3.17d) 


As these are (M+1 ) first order differential eauations (5 *17 )» (M+1 ) 
initial conditions are also needed. These are obtained from eqns(5 .16 ) 
by setting Fo = 0. 
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0 o (l»O)=1 (5.18a) 

9 1 (1,0) = 0 2 (1,O) = e M (l,0) - 0 (5.18b) 

and -which on substitution into eqn.(3*15) becomes 

6(1,0) = 1 (3.19) 


Values of factors e ^ and h^ in eqn. (3*17d), may be calculated using 
following two properties of the kernel 


M 0 - J [ f - e_] dFo - i (6 Jvl) 

o n=1 


( 3 . 20 ) 


M 1 = J [f -6 0 }Fo dFo = l (B/O 
o n=1 

Values of M 0 and M-j are given in Table 1 . 

One way to find out the values of and b^ is as follows: 


(3.21) 


assume b^ = U 


M 


( 3 . 22 ) 


Then 


J 


[fa - 3 3 dFo 


¥ S n Vm 

' nil »T + 7\T 


(3-23) 


Therefore, 


M-1 3 


I “ - )/ 3 


M M 0 n=1 Va 

and last differential equation (3*17d) changes to 


(3-24) 


d9 m (l,Fo) 

— 


Biref ( 1 » p °) }0 0> F< 3 M * M 9 M ( 1 * Po ) 


(3-25) 


- e 8 
MM 
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Chen, these simultaneous differential equation can he solved by any 
standard method. In the present study Runge-Kutta fourth order formula is 
used for solving these differential equations. 

3.3 HOMERIC ll Dr^FlCUlTIBS Ml THEIR REMEDIES 

3.3.1 SOTMA.TION OP INFINITE SERIES i 

Zorns 1 of integral' eon. ( 2 . 20 ) is a convergent infinite 
series. Series contains two factors, oscillating at different fre- 
quencies and their superimposed frequency is different from both. 

The factors are (-1 ) n and cos n ttx. 

For the numerical summation, above series has to be truncated 

at a suitable point such that terms after than point can be neglected 
in comparison to sum of the series. If this criteria is based on a 
comparison of single term, after the truncation point, with sum of the 
series, results will he in --error. The error arises because any single 
term may he small in comparison with the sum though the rest of the 
truncated terms may not he ne gligdble. The other reason is, individual 
terms may he laige enough, but the algebraic sum may be small because 
of oscillating nature of series. Hence the criteria fear truncation 
should he based on a comparison of sum of a complete cycle with sum 
of the series. 

3.3.2 CHOICE OP TIME STEP IN EQJT.(3.10)s 

A good plot between temperature and time is required for 
analysing a quench hardening. For a close to exact solution and 
smooth curve, temperature must be known after each small time interval 
(A t in eqn. 3 . 10 ). Choice of this time interval is very critical. 
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It depends upon the maximum Biot number. Too small a time interval 
will lead to a very large computation time while a laige time interval 
will lead to a divergence in the numerical scheme. 

Observing eqn.(3*10) for m=1 , equation gives a relation for 
9 (l ,Fo) as 

1- 0(l,Fo) “Bi ref h {9(1, T )}0 (l,T ) (At)^ (3.26) 


This relation can be reduced to 

C 1 


At 


r 


L Bi ref B{9 (1 ’ T ) }9 ^» T ) 


(5-27) 


where c^ is a constant of proportionality. Value of c.j is found by 
comparing results for different values of Biot number and time step, 

A t, with Heisler charts. This enables to make an optimum choice of At as 


A t 


= ChOI 

Bi ref S {9(1 , t ) > 9 (1,T ) 


( 3 . 28 ) 


3.4.1 EXIBRIMEHTAL SST-UPs 

Some experiment were also performed to verify the numerical 
results for cylinder and plate. A set up was made consisting of 
a furnace, a water reservoir, a recording unit and a transporting 
mechanism. The transporting mechanism was an assembly of a drum and 
a platform which could move on a pair of rails, (see photograph on page 31 )* 

The component can be held in any position with the help of a 
wire rope wraped on the drum. The recording unit was an Encardio-Rite 
four-channel recorder. It is made of the differential amplifiers 
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for amplifying the thermocouple output, the galvano meter for deflecting 
the needle and the chart drive motor with a variable speed gear set. 

The component -which made of mild steel was heated in the furnace 
at a constant temperature for a sufficiently long time to ensure the 
uniform temperature. The heated component was quenched by suspend ong 
it in the middle of water reservoir which was the quenching medium. 

The temperatures were recorded with the help of thermocouples embedded 
at center and at the surface. For ensuring the proper contact at the 
center, an arrangement was used as shown in Fig . 4 . A chart speed of 
1 mm/sec. was used for recording the temperature. 

The thermocouple used was Alume 1-Chrome 1 of gage 24 B&S insulated 
with Teflon. ISA code for these are KN and KP. These can measure a 
temperature upto 280°c with an error of +2$ and upto 1260°C with an 
error of +0.6$. 

34-2 EXEER MENTAL DIFFICULTIES s 

The main problem faced was to shield the thermocouple frcm water. 
For this purpose the Sauereisen insolute adhesive cement was used but 
this cement desolves in the water. Hence, a layer of Aral-dite was 
coated on the cement which can withstand temperatures of the order of 
200°C (used in the experiments). For the higher temperature, one has 
to look for a better coating. 

The results obtained by numerical method and experiments are 
compared in the next chapter. 






; 



CHAPTER 4 


KESUIffS AND DISCUSSION 

Solution obtained using numerical methods discussed in last 
chapter are presented here. First the numerical methods are 
compared with standard results available from Heisler's charts. 

This ccmparision will give some idea of the accuracy of the methods. 
Then these results are compared with the results available in 
literature and those obtained by experiments. 

In Table 2, the results of successive iteration method for 
the case of constant heat transfer coefficient (and various values 
of time stops) are compared with results of Heisler (1947) for 
various values of Biot number. It can be seen from Table 2 that 
for Biot nos. of 2,10,200, the optimum time step for which the 
agreement with Heisler's results are best are 0.01, 0. OOQ 4 ard 
0.000001 respectively. These results are in agreement with the 
dependence of time step on Biot number assumed in eqn. (3*27) and the 
value of C-j is seen to be 0.01. 

In Fig. 5, the results of successive iterative method for 
different Biot number are compared with Heisler's results. The 
results are in perfect agreement for low Biot number but becomes 
progressively worse as Biot number increases. It is also observed 
that error increases at large Fourier number. This can be inferred 
from equation (3»7)* The infinite series in equation ( 3 . 7 ) converges 
very slowly for small Fourier numbers. Any error in the evaluation of 
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Table 2 
Bi = 2 


Fo 

Heisler ’s 
results 


Time step 


O.OO5 

0.01 

0.05 

0.1 

0.55354 

0.55393 

0.55355 

O.54485 

0.2 

0.45759 

0.45792 

045770 

O.454II 

0.3 

0.35231 

0.35248 

0.35237 

O.55 066 

0.4 

O.2209I 

0.22090 

0.22085 

0.21969 

0.? 

O.19672 

0.19669 

0.19662 



Bi = 10 

Fo 

Heisler ’s 


Time step 



xesu i. us 

0.0001 

O.OOC4 

0.001 

O.O4 

0.25928 

0.26493 

0.26469 

0.5477 

0.08 

0.19467 

0.20232 

0.20214 

0.18543 

0.12 

0.16054 

0.17149 

O.17136 

0.14828 

0.16 

O.15837 

0.14973 

O.I496I 

0.12609 

0.20 

0.12429 

0.13653 

0.13642 

O.IIO96 


Bi ■ 200 

Fo 

Heisler *s 
results 


Time step 



0 . 5 x 10-6 

10" 6 

10“ 5 

0.0004 

0.26153 

0.27746 

0.27491 

0.25013 

0.0008 

O.19192 

0.21128 

0.20887 

0.17352 

0.0012 

O.I544I 

0.17578 

O.I7353 

0.13447 

0.0016 

0.12876 

0.15386 

O.I5I45 

O.IO38I 

0.0020 

O.II349 

0.13903 

0.13823 

- 
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the sum of the infinite series is amplified at large Biot numbers. 
Hence this method will lead to significant errors at high Biot numbers 
and low Fourier numbers. Note that at high Biot numbers, the cooling 
rate is high and hence the Fourier number of interest will necessarily 
be small. The error increases as time increases because of the 
accumulation of error. 


The separable kernel method is first validated for constant 
heat transfer coefficient by comparing it with results presented by 
Iieisler (1947 )• In Fig .6, it is observed that the agreement is again 
perfect for low Biot no. but results deviate more and more from the 
Heisler's results for high Biot number. It is also observed that 
separable kernel method initially overpredicts the temperature 
but after a particular time it tends to underpredict. The maximum error 
in temperature is found to be at low fourier numbers where it is about 

16$. But, if the error is based on the prediction time to reach a 

^ime-—... - Time,.. . 

SKM Heisler 


particular temperature (i.e. 


Time HMsler 


) , the error is 


considerably lower (of the order of 5 c /°)» 


3h Fig *7, the effect of the number of simultaneous differential 

equations (3*17) used on the accuracy of separable kernel method is 

seen. It is observed that as the number of equations are increased 

separable kernel method predicts the values approaching the Heisler's. 

It implies that more terms should be considered in the infinite series 

f f 'r high Biot number such that the last term of the kernel, 

e kp(-y 2 f Q ) becomes negligible, 
n 
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Ih Fig. 8, the results of successive iteration method for a plate 

for two different values of constant in eqn.(M*5)> are compared 

with the experimental observations. It Is found that results for 

C 8 £= 0.006 matches the experimental results better than that of 

O.OI 3 . This can also be seen in Fig -9 ^ J rhere the separable 

kernel method is compared with experimental results for two values 

of C sf The value of C depends strongly on the condition of the 

surface. Unfortunately, where the value of C - are given in literature 

S JL 

no information about surface roughness is given. In the present 
experiments, the mild steel surface had not been ground hence the 
surface was quite rough. Therefore, it seems justifiable to use 
the low value of C sf (which will result in higher boiling heat transfer 
coefficient, see eqn.(A4 *5 )• 

Successive iteration method may be called as exact method 
because ideally accuracy can be maintained by a proper choice of time 
step and more number of iterations. But, if it is not possible to 
avoid round-off error (because of limit' * on the number of significant 
digits that can be stored on the computer) then this method is not 
useful since the error tends to accumulate with time . On the other 
hand, in the separable kernel method the error in the previous step 
damps out rapidly. This is because the error in Runge-Kutta method, 
which is used in solution of simultaneous differential equations ( 3 . 17 ), 
decreases if the slope of the temperature - time curve decreases, 

(see Rajaraman (1974)). In addition the accuracy can be increased to 
any desired degree by simply increasing the number of equations. This 
is illustrated in Fig. 10. 
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In Pig.11, experimental results for the temperature at the plate 
centre are compared with that of separable kernel method. It can be 
seen that numerical scheme shows a continuous fall in temperature 
from the very beginning, while experimental results show no fall for 
sometime. The variation of temperature at the centre is slower in 
the beginning and it looks as though the recording instrument is 
unable to sense it. At large f ourier number error is much higher.* 

Same trend is observed for cylinder surface and centre 
also (Pigs. 12,13). 

In the above presentation film boiling could not be 
illustrated. This was because of laboratory limitations, it 
was not possible to reach high temperature needed to achieve film 
boiling. However, experimental results for film boiling are 
available from Bofors Handbook (Thelning 1975) and these compared 
with the numerical solution in Pigs. 14,15. Though the deviation is 
very large, trend is again same as in the previous illustrations. 
Temperature predicted is higher than the experimental results 
available . The reason may be an improper choice of correlations. 
Unfortunately, the conditions are not mentioned in the literature, also 
he mentioned that the cooling capacity of water is increased very 
considerably by adding 10$ of common salt or soda to the water. 

With such a information, the over-prediction of temperature was 
expected. 

* One reason for higher error may be thought as end losses due 
to limited size of specimen. Bat the error was, calculated by 
so called Newman method, very small because of the low Fourier number 


for the ends 
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CONCLUSION : The present work is not merely of theortical interest. 

In the industries , often we. face an inverse problem, where one desires a 
conponent of specific alloy with specific qualities and suitable rate of 
cooling has to be found cut. The present computing techniques will judge 
whether the choice of cooling system is proper. If it is not so, one 
may vary the cooling rate by using another quenchent or adding seme 
chemicals to the same quenchent* Cooling rate may also be increased by 
mechanical vibration to the component in film boiling regim. This inter- 
mittently breaks the film layer. 

We may control the depth of hardness also by choosing the proper 
cooling rate* From the figs. 10 and 11, we can see that, the cooling rate 
at centre is faster than at the surface* Hence, by increasing the cooling 
rate we can decrease the depth of hardness and vice-versa* 

Keeping the dignosis of error in two methods in view, we 
may use the two methods simultaneously for more accurate results. 

The successive iteration method may be used for low fourier number, where 
one needs more accuracy, without much increasing computation time* 

But at the later stage, when successive iteration method begins to 
deviate due to accumulation of error, we may switch over to separable 


kernel method* 
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AFBBHDIX I 


The following expressions are known as Theta functions. 

m 

2 


0,(y ) « 8,( p/x ) = l exp(iir t n ) exp(2npi) 

J J n=-QD 


(A1 .1 ) 


and 

V u ) 


V ' 2 

1+2 lexpCi^Th ) cos2np 
n=1 


6/(y/ T )= I (-1 ) n expCi^x n^) ■ exp(2nPi) (H.2) 
n=“co 


- 1+2 I H) n e 1 7r1h c0s2n y 
n=1 


Tiiere exp(i TTt) = q 


and 


1 I < 


(A1 .3 ) 
(A1 .4 ) 


Hence, singular kernel of eqn.(2.20) is also a ^eta function 


1+2 l (-l) cos nir x exrf -n tt (fo-t)} 
n=1 


= 9 4 {“^ / i HFo-t ) } 


6 5 / i TT (po- T )} 


(A1 .5) 


Some relations are given "by Gradshteyn and Eyzhik (1965 )* few of those 


are mentioned below. 


e 4 (y + w) - © 4 (p ) = 0 4 (~u ) 

(ai.6) 

0 5 (V + w) =0 5 (u ) = 0 3 (-V ) 

(A1 *7) 

® 4 ( P+ ^ ) =6 ? ( W) 

(A1.8) 

w+ jr) = 0 4 ( v) 

(A1.9) 
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APBSM5IX II 

INVERSE LAHACE TPJUTSFOMATIQN BY SERIES EXPASFS I® 

In this method of Inverse Iaplace transformat ion.', F(s), 

the image function of f ( t ) is expanded in a series, 

00 

P(s) = I F (s) (A2.1) 

n==c 

in which the terms F^(s ) are image iscnfunctions of f (t) * Then solution 
is obtained by transforming this series term by term, Doetsch ( 1 97 ^ ) 

CO 

f(t)« I f n (t) (A2 .2 ) 

n=o 

This cannot always be correct because it amounts to the inter- 
change of an infinite summation and an integration. But certain types 
of series can be transformed without any error. Che of those is described 
below • 

If F(s) is a rational function 

r l(s) 

(A2.3) 

in which x^(a) is. of lower degree than the denominator ^(s), then 

F(s) can be expanded In partial fractions. If denominator has the 

s Imple zeros at ou j ct ^ »ot , ...,a , ... 

* ^ o n 


r 2 (s) = (s- 

a-j )(s-a 2 )(s-a ^)...(s-a n ). .. 

(A2.4) 

IT 

F(s) « l 
n- 1 

r F V i 

(42.5) 



r l( a n > 
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where 


!(° J 


r 2 


are residues at poles a * 

n 


How, transf ormation is 


made term by term and the result is 


f(t) = 


H 


n 


r 1 < a n^ 


exp( a t) 


'" n 


(A2 .6) 


le place inverse transform of function 

I_(r /s) 

P(b) = -Si i 

✓s X,(/s) 


(12.7) 


can be obtained using above mentioned series expansion method. The 
function P(s) of eqn.(2.7) has S explicitly in its denominator, 
therefore as s -»■<»» P(s)-*- 0. One pole of function F(s) is at zero 
and other infinite poles are given by 


I 1 (/s)=-J- J^(i /s) = 0 


(A2.8) 


or 



(A2.9) 


where, U n is root of equation 


Vuj-O 

Then residues at poles are given as 

*l( <^) ^ s ) 


i(a ) 
>v n' 




(A2.10) 


Taking derivative 

r 1 (“n> 


E^(a n ) 


x o(^ ’ / s) 


1 X,( /s) +§■ IJ( /s)} 


2/s 


s *“*4 


( 12 . 12 ) 
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Now, since 


I Q (r / s) 

s+ -v 2 n 

= I 0 (ir/ s) 

s-> -y n 

- J o< »n 

(12.15) 

If(/s) 

n 

= J^*(i 4) 

s ^y n 

J i< V 

( 12 . 14 ) 

and Xj ( / s ) 


= 0 






Above mentioned eqn.(A2.12) for residues reduces to 


r l( a J 

r ^(a n ) 


2J 0 ( ^ r > 
»n> 


(12.15) 


Using recurrence relation we obtain. 

J 1 ^«) 

J AV = J o<“n>- — 


(12.16) 


n 


Since, according to eqn. (12 .10), J^(n ) = 0, eqn. (12. 15) reduces to 

r lU„) 2J 0 U„ r ) 


and 


r 2< %) 


r-j (o) 


r 2 '(o) 


J o< V 


2J Q (ix /s) 
J 0 (i ^s) 


(12.17) 


(12.18) 


s=o 


Hence, liaplace inverse of function f(s), eqn.(A2.7) can be written in 

the following form with the help of eqns. (A2.6)j(A2.17 ) and (12.18). 
(o) » ^(a ) 

f ( t ) = — + l — r . exp(a n t) 


vj(o) n-1 a a ) 


v J o ^ n / ,,2 ^ 

2+2, 2 ——7 r~ exp(~u t; 

n=1 Vun j n 


^ X'12.1 

GH«m; -ism 


Aoc. No. i“'\ 


54916 
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APBSND3X III 


Eie following formula is known as Poisson summation formula 
Giver (1974), Askay (1975 5‘, 

CO CO CO 00 

£ 1 f(n) = J f(n)dn + 2 £ J f(n) cos2^s n d n (A3«0 
n=o o s=1 o 


Prime, as a superscript of summation sign, means value of f(n) is 
half of it, at n = o. 


Singular kernel of eqns. (2.20) is given as 

CO 03 

i t* ii 2 2 

l f(n) = 1+2 I (-1 ) cos n^x exp{~n tt (fo-t)} 
n=1 n=1 


OOf 

= 2 £ (~1 ) n cos n^x exp(-n^TT^(Fo- T )} (A3.2) 

n=o 


On the application of Poisson summation formula 


CO 00 p p 

£ 1 f(n) = 2 [ J (-1 y 1 cos n-n-rx exp{-n ^ (I'd- t) }dn 
n=o o 


xi 2 2 

J (-1 ) cosn ttx exp {-n ir (Po- t) }cos2-nrsn dn] 
o 

(A3. 3) 


+2 


Ji 


Gn evaluating the integrals 

1 


I 1 f(n) = 2 I 
n= 1 


l(. 


4/fr(Fo- t)) 
1 


, . ( x-l f. , 

{ exp{ A + exp { 


4(Fo~ t ) 


4(Fo~ t ) 


»} 


+2 


I < 

s=1 4 /{ , h(Po-t)} 


(x+2s-l) 2 (l+2s-x) 2 ~ 
{ exp( — _)+exp( — _) }} 


4(Fo-x) 


4(Fo-t) 

(A3 .4) 
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and then rearranging, we get 

I’ f ( n ) = — ; A exp{ -(x-1) 2 /4(Fo- x ) 

n=1 /{ tt ( Fo- x ) } 


+ l exp { -(x+2s-1 ) 2 /4 (Po- t )} + exp { -(l+2s-x) 2 /4(Fo- 

s=1 


(A3. 5) 


For the surface, this kernel reduces to 

CO 

I CO 

I f(n) = [1+2 l exp { -s / (Fo- t ) } ] {ir(Po- t ) > 2 
n=1 s=1 

Replacing s by n, eon. (A3 .6) can be rewritten as 

CO CO 

1 f(n) = [ 1+2 I exp{ -n/(Fo- t)}] {it (Fo- t) 
n=1 n=l ' ' 


(A3 .6) 


T )>} 


(A3. 7) 
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APPENDIX IF 

Computer programs for successive iteration method and separable 
kernel method were developed for numerical solution of the nonlinear 
singular Voltsrra integral eqjuation. The program is generalised for 
three simple geometries, i.e. plate, cylinder and sphere distinguished 
by a variable name IGMTRY. 

The computer program is written in Fortran IF. Fortran listing 
of this program is given in Appendix F. 

1. Description of Programs 

MAIN program reads various inputs , calls subroutines for surface 
and centre temperature (step by step) and prints the output data. 

The two methods differ in only the way they calculate surface 
temperature. For separable kernel method subroutine DJiRIF is called from 
subroutine SURSKM while, for successive iteration method subroutine 
GITA is called from subroutine SURSM. 

Subroutine Name Description 

Xt reads all the input, calculates the time 
step, calls subroutines EPS, SbRSKf.^/ STIRS HI and 
CENTRE and prints the output. 

STJRSKM Used for calculating surface temperature. It 

(used in separable solves the set of simultaneous differential 

kernel method) . 

equations (5.17) by fourth order Runge-Kutta 

method* It calls subroutine DERIVE* 
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DEBIT 
SUES 3M 

(used in successive 
iteration method) 

GITA 

EPS 

CENTRE 

BESJ 

KNOW 

FLUX 

Program Symbol 

ALPHA 

BETA 

BETAO 

BETAK 

UMAX 

CPF 

CPS 

CPV 


It calculates the derivatives. 

Used for calculating surface temperature 
by Poission summation formula. It calls 
subroutine KNOWN and GITA. 

It calculates the last integral in eqn.( 
by the Modern Gaussian integration formula. 

It calculates the eigenvalues and other 
constants. 

It calculates the temperature at the centre 
and for that calls subroutine KNOW. 

It calculates the value of Bessel's function 
of zeroth order. 

It evaluates the integral by Simpson's 
rule . 

It calculates the outgoing flux corresponding 
to a nondimensional temperature. 

Definition 
Thermal Diffusivity. 

Coefficient of thermal expansion 

® o 
0 

k 

Max. Biot number. 


Thermal capacity of fluid 
Thermal capacity of solid 
Thermal capacity of vapour 



Constant of Nucleate boiling Heat transfer 
correlation . 


D- 

Thickness of solid 

DELT 

Small time step 

EIGEN 

Eigenvalue 

EMSVTf 

Emissivity of metal 

EPSN 

£„ r } a constant ' 

FLITS 

Outgoing flux 

FONO 

Fourier number 

G 

Gravitational constant 

GC 

Conversion factor 

H 

Film Heat transfer coefficient 

HB 

Boiling component in the film boiling 

UMTET 

Geometry of solid 

KF 

Thermal conductivity of fluid 

KS 

Thermal conductivity of solid 

EV 

Thermal conductivity of vapour 

L 

Height of solid 

MCJF 

Viscosity of fluid 

MJV 

Viscosity of vapour 

N 

Number of equations used. 

PI 

ir 

PR 

Erandtl number 

RHOF 

Mass density of fluid 

RHOS 

Mass density of solid 



RHOV 

sk-mai 

SIGMA2 

T 

TP 

THETA 

TBSTAC 

THETAI 

THETAI/ 


THETAS 


TH3TAT 


TIHTL 

TLEIDH 

TSAT 

TTRAHS 

X 

Z 


Mass density of vapour 
Surface tension 
Steffen Boltzman cons t&nt 
Temperature component 
Temperature of fluid 
Hondimensional temperature. 

Hond-imensional temperature at centre 
Hondimensional temperature corresponding to 
initial temperature 

Hondimensional temperature corresponding to 
Leidenfrost temperature 

Hondimensional temperature corresponding to 
saturation temperature 

Hondimensional temperature corresponding to 

transition temperature 

initial temperature 

Leidenfrost temperature 

Saturation temperature 

Transition temperature 

Hondimensional value of coordinate 

Value of integral. 
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2. Correlations used in the numerical solution; 
a. Film Boilings 

In the film boiling regim heat transfer by radiation dominates 
over heat transfer by convection. 


Boiling heat transfer coefficient h b is given as 


h b = O.943 


! K 5 P ( 


i v V' r 


p . - p ,)s x ' 


! 1 Ur 


■where 


X* = h- (1 + 

fg 


0.4 (T - tJ C, 


f' u pv 


h 


fg 


(M-1) 


(M.2) 


and total heat transfer coefficient h is given by 

h b 1/5 

h -h(-> +h r 


(M- 3) 


here h r is heat transfer coefficient by radiation given as 

(t! - 4) 


where 


h r = 0e 


P v 

Pf 


T 


T, 


s ~f 

Thermal conductivity of vapour 
Vapour density 
Fluid density 
g = Gravitational constant 
L = Height of the surface 
= Viscosity of the vapour 
C = Thermal capacity of vapour 

if » 

h^g= Latent heat of vapourisation 
a - Steffen Boltzman constant 
e = Emisf.sivity 
T = Surface temperature 
Tf- Fluid temperature. 
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b.- nucleate Boilings 

Correlation for nucleate boiling heat transfer is 


h f g P r 1Tr 


= c 


j h(T s -T f ) 


sf 


./■ 


s c cr 


! 


f n fg ' v/ P f ~ P v'L 


0.53 


"Pf 


J sf 


There = Thermal capacity of fluid 

A constant 

U = Viscosity of fluid 

w = Conversion factor 

°o 

a = Surface tension 


c. itee Convections 

Che correlation for free convection is 


HU - O.555 (GrPr) 


tjhere Hu = Suss e It number 
Gr = Grashof number 


(M -5 ) 


(A4 .6 ) 


Rr = Prandtl number 
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3* Input data for numerical solutions 

a. Vapour properties (at 100°c)s 

Density 

Thermal conductivity 
Thermal capacity 
Viscosity 

latent heat of vaporisation 

b. Fluid properties (at I00°c)s 
Density 

Thermal conductivity 
Thermal capacity 
Viscosity- 
Surface tension 


O .593 %/m 5 
0.0249 j/m-°K-sec. 
2.028 Kj/kg~°K 
0.0121 mil Sec/m 2 
2256.94 Kj/% 

959 %/m 5 
O .669 j/n~°K-Sec 
4.2:6 Kj/%-°K 
0.278 mN Sec /m 
58.9 mIT/m 


c. Fluid properties (at 70°c)s 

(considered at the mean temperature in free convection regim) 
Density = 990 %/m 5 

O 

Viscosity = O .57 mlT Sec/m 

Coefficient of Thermal expansion = 0.4x10 5 l/°K 

Prandtl number = 3.6 


d. Material properties* 
Dens ity 

Thermal conductivity 
Thermal Capacity 
Bnissivity 


= 7760 %/m 5 
= 36 J/m-°K-Sec 
= O .486 Kj/%-°K 
= 0.9 


e. Specimen data* 

Plate 

cylinder 

Fluid temperature 

Initial temperature for 
experiments 

Initial temperature for the 

results in Bofor's Handbook 


25 x 25 x 3.7 cms. 

20 cms. high, 5*0 eras. diameter 
= 26 °C 

= 184% 

= 900%. 
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000 ooo ooooo ooo 


MAI N PROGRAM 


REAL KV,KF, KS, L,MUV, W LF 

f* OH A k‘M Oi ?i i T » r i r* a t % 


COM MO I’ 
CCMMCf 


K1-, Kb, L , MUV » V L r 

(5) , FOND l 1O0D }, THETA ( 1 COO) » THE TAC , K , DE LT 
IG C H(50 ), EPSN, BETAC, B2TAK 
W/X f 16 ) ,W { 16 ) 


COMMON E I GEN { 50 3, EPSN, EETAC, BETAK 

COM MON/ XW/X ( 16 ) ,W { 16 ) 

CCXX0N/CQ'!ST/C1,C2,C3,CA,C5,C6,C7,B 

CCXMO.'J/TE'-'X/TINTL , T F , T FETA L , THETA T» THETAS 

h PCR NO. IF ECHOS. IN SKM.NO. OF GAUSSIAN ORDER IN SIM 

METHOD =1 SUCCESSIVE ITERATION MFTHHE (SIM) 


N PCR NO. 
METHOD =1 
METHOD =2 
READ XO , 
I GM TRY -2 
IGMTRY =1 


s fvlfvAIUJiUifU f U U f U f ■f LJ 

X/TIN'TL , T F, THETA L, THETA T, THETAS 
IF ECHOS. IN SKM.NO. OF GAUSSIAN I 
SUCCESSIVE ITE R AT I 0. . METHOD (SIM) 
SEPARABLE KERNEL METHOD { SKM ) 

!, METHOD 

INFINITE CYLINDER . 

INFINITE PLATE- 


INPUT DATA 


4 CONTINUE 
PRINT 500 

5 READ ICO, IGMTRY 

I F ( IGMTRY-i ) 70,6,7 

6 PRINT 601 
GC TC 3 

7 PRINT 602 

8 CONTINUE 

READ 400, T I NT L ,TL E I CN, TTR ANS , TSAT , TF 

READ 400, R HCV , KV , MUV, CPV, HFG 

READ .400, RHOF, KF, MUF, CPF, BETA, SIGMA 1 

READ 400 » RHOS » KS » CPS 

READ 400, G,GC,CSF, S IGMA2, EM SVTY 

READ 400, L,D 

EVALUATION OF CONSTANTS 

CALL EPS ( IGMTRY, N) 

THE TA I = (T I NTL-TF ) / ( T I NTL-TF ) 

THE TAL= (TLE I CN-T F ) / { T I NTL-TF ) 

THETAT=(TTRANS-TF)/{TIML-TF) 

THETAS=(TSAT-TF)/{ T I NTL-TF ) 

ALPHA=KS/{RKCS*CPS)' 

PR= CPF*MUF/KF 
Cl= HFG 

C2=0.4*CPV* (TI NTL-TF ) 

C3= 0.943* (KV**3*RH0V*( RHOF-RHOV) *G/ t L*MUV* < TIN TL-TF ) ) )**.25 
C4»S,IGMA2*EMSVTY 

C5= UCPF/(CSF*HFG*PR**1.7) )**3 )*l TI NTL-TF ) **2*MUF*HFG*SQRT (G* ( RHOF 
1 -RH0V)/(GC*SIGMA1) ) 

READ 400, PR, RHOF, BETA, MUF 

C6= .021*(P n .*RHCF**2*BETA*G*ITINTL-TF) /MUF**2 )**0.4*L**0. 2 
C7= ( D/2. ) /K S 
PRINT 700 

INITIATING THE CALCULATION 


TT= 1. 

J=2 

GTT=FLUX { J,TT) 
BMAX=B 

DEL TIN = .001 /BMAX 
DEL T= CELT IN 
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FC ■] C ( 1 ) = DELT 
TFET A ( 1) = 1. 

TF FT AC = 1 • 

T ( 1 ) = 1 » 

CC ID 1 = 2. • 
l'< Til )= '. 

K = '> 

2 ' CC ,T INL V 
- K=K+1 

IF{ vPTI-OC.rC .2 ) GC TC 21 
CALL SIRSPMM 
T5MP=TFFTA( K ) 

CC TO 31 

21 CALL SLFSK-HM 
T£MP=;;. 

CC 30 1 = 1, M 

30 TEM P=T FMP+ T { I ) 

THETA ( KJ=TEMP 

31 CCMTINUF 
CALL CENTRA 

40 CONTINUE 

PRINT 200, FCNC(K),TFFTA(K),THETAC 

5 ) C C\‘ TINU C " 

IF(FONOtK). GE. 5 . ,CR .TEMP.L E. . i) GO TO 60 
GTT=FLUX{TFMP ) 

CELT=CELTIN*{BNAX/3 )**2 . ■ 

FONO ( K ) = FH'3 0 ( K~ 1 ) +DELT 
GC TO 2 i 
6 j CONTINUE 
GC TO 5 
70 CONTINUE 
STOP 

10 ) FORMAT f 2 >I4) 

2C) FORMAT {1 OX, E 12 . 5 , 5X , FF • 5, 5X, F 8. 5 ) 

300 FORMAT (//F12. 7, 14//) 

400 FORMAT (6E12.5 ) 

5 GO FORMAT (1H,*, GENERAL PROGRAMME FOR SURFACE TEMP. FOR ANY GEOMETRY 
1Y BY SEPARABLE KERNEL METHOD*) 

601 FORMAT (//1DX»*GE0METRY — INFINITE PLATE*/) 

602 FORMAT (/ /ID X,*GFOMETRY — INFINITE CYLINDER*/) 

700 FORMAT ( 16X, *FONQ* , 7X , *TEM P S* , 5X , *TE MPC*/ / ) 

END 
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20 


SURSKM (N ) 

COMMON T ( 5 ) , FONO ( 1000 ) , THETA ( 1000 ), THETAC, K , DELT 
COMMON EIGEN (5C ), EPSN, 8ETAC, BETAK 
DIMENSION S ( 5) ,S1(-5),S2( 5) ,S3( 5) , S4( 5) 

DO 10 1=1, N 
S C I ) = J. 


CALL DERIVf N,0»,S,S1) 

CALL DERIVf N,. 5, SI, S2) 

CALL DERIVfN,.5,S2,S3) 

CALL DERIVf N»1.,S3,S4) 

DO 20 T — 1 f\| 

T ( I )=T ( I ) +OELT * { S If I ) + 2.*S2( I ) + 2.*S3f I ) + S4(I))/6. 


RETURN 

END 



S U 3 C HUT I ;•/ * C C R IV { N» C, S , SS ) 

CCM MCN T { 5 5 ,FONC( 1000) , THETA { 1000), THETAC,K,DELT 
C CM VON r I G ' ' J ( 5 C ) , EPSN , BETAC,8ETAK 
DIMENSION c { 5 ) , S S ( 5 ) 

TT=\'. 

DC 10 1=1, ' 

10 TT=TT+T( I) + S( I }*DELT*C 
GTT = FLUX (TT ) 

X=FCfJO(K )+OELT*C 
Y=T ( I ) +S ( I )*C=LT*C 
DC 60 1=1," 

IF(I-l) 60,3'', 21 
2 • I F f I— M ) 40, 5^, 60 

3 ’ SS( I ) =-R ETA i*C-TT 

GC TO 60 

4 * SS( I )=~6ST4K*GTT-EIGEI\ ( 1-1 )**2*Y 

GC TC 60 

5) S$( I ) = -EPS : ’ * p ET AK^GTT— £ I GE M I — 1 ) **2 *Y 
6i CCNTINIF 
RETURN 
END 


SUBROUTINE -PS ( I G HTRY , N ) 

COMMON T ( 5) , FOND { 1 10 0 ), THE TA ( 1000 ) , THE TAC , K , DE LT 

CCMPCN E I O r ! ■! ( 50 } » EPSN , BET AO, BETAK 

N1=N-1 

N2- N— 2 

GO TO (100,200,300), IGHTRY 
100 BET A0= 1. 

BET AK = 2* 

P 1= 3 • 1415^^ 7 
CO 110 1=1, Ml 
11' EIGEN{ I)=PI*FLCAT( I ) 

EPS N=1 • /3» 

GO TO 400 
200 BET A0=2. 

BET AK = 2. 

EIGENC 1)=3. 8317 
EIGEN ( 2) =7. 0156 
EIGEN (. 3) = 1D*1725 
E.IGEN(4) = 13*3227 
DC 210 1=5, N1 

210 EIGEN { I ) =FLQAT ( I)*PI+PI/4. 

EPSN= « 25 
GC TO 400 
300 CONTINUE 
400 CONTINUE 

CC 420 1=1, N 2 

420 EPS N=EPSN-B ETAK/E I GEN( I)** 2 
EPSN=EPSN*EIGEN(N1)**2/BETAK 
RETURN 
END 
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SIHROUTK- StJRSIM(H) 

COMMON T ( 5 J , F 0 N 0 ( 10 K‘ } , THETA ( 1000 ), THE TAC , K , DE LT 
COMMON E IGFN (50 ) , E°SK» BETAC, BETAK 
CC !«CL/XW/T { 16 ) ,W( 16 ) 

CI'^ENS I<'N »CIS»’M( 16 5 
-EAL I NT 6 U ’ 

TH' : T AH = 1 « 

IFf K.ME.l) "F~7AN-TFET A { K—i ) 

CO 12.' i = i,;; ■ 

s= *. 

M=i 

110 CONTINUE 

PCI SUM { I )=S + FXP{-PLOAT(M)**?/<OELT*X{ I 1**2) ) 

IFf ABSfPCISUMf I)-S).LT..C01) GO TO 115 
S = ° C I S U M { I ) 

M = M + 1 
GO TO 11'- 
115 CCITIK'tJF 
120 CCOTIUUF 

C TAKING FIRST APPROXIMATION OF THETA ( K ) =THETAM 

THE T A ( K ) =TH ~T AN 
CALL KNOWN (Z) 

M=1 

3 CO C O’ 1 T I f ! U E 

CALL G IT A ( '■ I , I NT FUN, PC I SUM ) 

S=Z - I NT FUN 

IFf ABS (TFE T A(K )-S J.LT. f . C C C<U*TH ETA { K } ) GO TO 400 
THETA (K)=S 
M = M + 1 
GC TO 300 
400 COUNT INIJE 
RET LRU 
END 
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SUBROUTINE GITAfN, INTFU’N, PGISUM) 

COMMON T (5) , FONOf 1000 ) tTHETAI 1000 ), THE TAC , K , DE LT 
COM MON FI GEN (50 ), EPS N, BET AO, BETAK 
COMMON/XW/X { 16 ) , W ( 16 ) 

REAL INTFUN 
J = 1 

PI=2.*ARSrif 1. ) 


G=F LUX ( J ,TEMP ) 

H=H + W { I)*(l .+2.*P0ISUMU ) )*Q 
CONTINUE 

INTFUN=2.*$QRT{ DELT/PI )*H 

RETURN 

END 


SUBROUTINE CENTRE 

CCMMON Tf 5) , FONOf 1000 ), THETA f 1000 ) » THETAC,K ,DE LT 
COMMON EIGFNI50 ), EPSN , BETAC, BETAK 
CALL KNOWN f Z ) 

RETURN 

END 



S 15 POUT I N T KNOUN (Z) 

CCM,vON 7(*) , FOLD { 1030 ) f THETA { 1000 ) , THE TAG, K 

CCMN'OO F? I G~ (50 ) , EPSN, BETAC, 8=TAK 
M =K 


X = ' ! , 

J=2 

p I=:. • * A P. 0 T 1 . { 1 • ■ 

z=o . 

I F ( W. LE. 1 } GO TO 60 

,W 1=^-1 

cc E*t r=i,'u 

Z 1= 1. 

I FC 1-1 ) 60, 10, 20 
10 TAU = 0.5*F0'!Gm 
GC TC 21 

20 T AU = FOND ! I ) ~ ! 'm 5*( FONC( I )-FOHO( 1-1) I 

21 CC‘TTIfJLF 
K = 1 

Z Z= -1 . 

20 CQJ T I HUE 
El= EIGEN (K) 

I F f I GM TRY-2 ) 3 ,31,22 

30 CONTINUE 

Z 2= Z 1 + 2 « * Z Z * EX P ( - E I * 2 * { FCNO { M ) - T AU ) ) 

GC TO 32 

31 CALL BESJC"I,n,3J,,091,IER) 

12=11 + 2 , *(1 ,/B J)*FXP(-EI**2*( FOMO(M J-TAU) ) 

32 CENTIME 

I F ( ABS (Z2-Z .1 ).LT. ( .01*Z2> ) GO TO 40 
Zl= Z2 
K=K + 1 
ZZ=~ZZ 
GC TO 20 
4n CENTIME 

I F ( 1-1 } 60,41,42 
T EM P= ( l.+THETA (1) ) / 2 . 

GC TO 43 

TEMP=( THETA ( I- 1 ) +THETA ( 1 ) ) /2. 

CONTINUE 
G=F LUX ( J , TE MP ) 

Z*Z+Z2*G*DFLT 
CONTINUE 
CONTINUE 
Z=1.0-Z 
RETURN 
END 


41 

42 

43 


50 

60 


,DELT 



2 > 


SUBROUTINE URS J (X«N»BJ tO* IER) 
SUO*OUTINE FCP- ZEROTH CRDEP C .LY 
IFC A.NE.O) GC TO 20 
X2= ( X/2* ) Y* 2 
B J = I # 

X3=l. 

K. =1 

X 3= — X 3 # X ? / F L C A T f x )=’--* 2 
X4= X3 
K = K + i 

X3=-X3*X2/FLCA7(K)#*? 

X4= X4 + X3 
BJ=EJ + X4 

IF{ ABS (X4). LE.AESt C*BJ ) ) GC TO 20 

K=K + 1 

GC TG 10 

CONTINUE 

RETURN 

END 


FUNCTION R, UX(J,TT) 

CCMW0N/C0'iST/Cl,C2,C3,C4,CStCfc,C7,B 

CCN IXCN/T EMX / T I N T L , T F , r f-ET AL , THET AT, THETA S 

REAL L EM CA 

GC TO (1 , 2 ! , J 

TT= 1. + (THETA (K )-l. )*{ 1 .-XU J**2) 

IF( K.l'iE.l) TT= T FETA{ K-l )+TFFTA(K )-THETA ( K-l ) )* {1. 
CONTINUE 

IF( TT.LE.THFTAL) GO TC 10 
LEMCA=C1+C2 *TT 
H B= C3 * { L E MD A / T T ) # *0 • 2 5 

HR=C4* ( (TT* (TIML-TF ) + TF + 273. ) **4- ( TF + 273. )**4) 
H=HE+HP 

H=HB*{ HB/H) **( 1./3. )+HR 


B=C7*H 

FLUX=8*TT 

RETURN 

10 IF( TT. LF.THETAT ) GO TC 20 
GC TC 20 
15 CONTINUE 
RETURN 

20 IFITT.LE. THETAS) GO TC 30 
H=C5*TT**2 
B=C7*H 
FLU X- B.*T T 
RETURN 
30 CCNTINUE 

H»C6*TT**ft.4 . ' ■ , 

8*C7*H. 

FLUX=E*TT . 

RETURN 

END 


-X (I )**2) 



